////////////////////////////////////////////////////////////////////////////////
// Options
////////////////////////////////////////////////////////////////////////////////

// To log all computation in a log file
#define LOG_COMPUTATION

////////////////////////////////////////////////////////////////////////////////
// Includes
////////////////////////////////////////////////////////////////////////////////

#include "Plugin.hpp"
#include "macros.h"

// Conversion functions
#include "Convert.hpp"

// Toolit
#include "TKMacros.h"

// IO stream
#include <iostream>

#ifdef LOG_COMPUTATION
	// STL classes
	#include <vector>
	#include <utility>
	#include <string>
	#include <fstream>
#endif

// Namespaces
using namespace SC;
using namespace std;
using namespace Toolkit;

////////////////////////////////////////////////////////////////////////////////
// Experiment data collection
////////////////////////////////////////////////////////////////////////////////

#ifdef LOG_COMPUTATION
	//static unsigned int nb_computations;
	static float model_time;
	static vector< pair<float, unsigned int> > quantity_A, quantity_B;
	static unsigned int bound;
	static string file_name;
#endif

////////////////////////////////////////////////////////////////////////////////
// Plugin init and exit
////////////////////////////////////////////////////////////////////////////////

void InitPlugin (const string & prog_name) {
#ifdef LOG_COMPUTATION
	//nb_computations = 0;
	model_time = 0.f;
	bound = 0;
	quantity_A.push_back(pair<float, unsigned int>(0.f, 0));
	quantity_B.push_back(pair<float, unsigned int>(0.f, 0));
	file_name = prog_name;
#endif
}

void ClosePlugin () {
#ifdef LOG_COMPUTATION
	quantity_A.push_back(pair<float, unsigned int>(model_time, quantity_A.back().second));
	quantity_B.push_back(pair<float, unsigned int>(model_time, quantity_B.back().second));

	// Open log stream
	string log_name = file_name + "_log.m";
	ofstream log_stream(log_name.c_str(), ios_base::out|ios_base::app);
	if (!log_stream.is_open()) {
		cerr << "Cannot open " << log_name << endl;
		return;
	}
	// Write data to file
	time_t tt = time(NULL);
	char buffer[32];
	strftime(buffer, 32, "%y%m%d_%H%M%S", localtime(&tt));
	log_stream << "bgn_" << buffer << ".a.x = [";
	for (unsigned int i=0; i<quantity_A.size(); ++i)
		log_stream << quantity_A[i].first << " ";
	log_stream << "];" << endl;
	log_stream << "bgn_" << buffer << ".a.y = [";
	for (unsigned int i=0; i<quantity_A.size(); ++i)
		log_stream << quantity_A[i].second << " ";
	log_stream << "];" << endl;
	log_stream << "bgn_" << buffer << ".b.x = [";
	for (unsigned int i=0; i<quantity_B.size(); ++i)
		log_stream << quantity_B[i].first << " ";
	log_stream << "];" << endl;
	log_stream << "bgn_" << buffer << ".b.y = [";
	for (unsigned int i=0; i<quantity_B.size(); ++i)
		log_stream << quantity_B[i].second << " ";
	log_stream << "];" << endl;
	// Close stream
	if (log_stream.is_open()) log_stream.close();

	// Open res stream
	string res_name = file_name + "_res.dat";
	ofstream res_stream(res_name.c_str(), ios_base::out|ios_base::app);
	if (!res_stream.is_open()) {
		cerr << "Cannot open " << res_name << endl;
		return;
	}
	// Write data to file
	if (quantity_A.back().second > quantity_B.back().second)
		res_stream << "1\t0" << endl;
	else
		res_stream << "0\t1" << endl;
	// Close stream
	if (res_stream.is_open()) res_stream.close();
#endif
}

////////////////////////////////////////////////////////////////////////////////
// SC Callbacks
////////////////////////////////////////////////////////////////////////////////

#define A_NORMAL_CREATE_RATE	0.20f

#define B_NORMAL_CREATE_RATE	0.37f
#define B_INHIBIT_CREATE_RATE	0.027f

#define AB_BIND_RATE			0.72f //* 0.42f
#define AB_DEGRADE_RATE			0.53f

#define A_DEGRADE_RATE			0.002f //0.0085f
#define B_DEGRADE_RATE			0.002f //0.034f

#define bA_BIND_RATE			0.19f
#define bA_UNBIND_RATE			0.42f //0.12f

bool CreateProtein (SCSystem & gene, const SCSystem &, SCSystem & protein) {
	// Get inhibition state
	bool inhibited = (gene.GetSchema1(INHIBITED_OFF, INHIBITED_LGH) == "1"),
		 type_a = (gene.GetKernel(K_A_OFF, K_A_LGH) == "0");
	// Get rate
	float rate = type_a ? A_NORMAL_CREATE_RATE : ( (!inhibited) ? B_NORMAL_CREATE_RATE : B_INHIBIT_CREATE_RATE);
	// Generate protein based on a rate
	if (frand() <= rate) {
		unsigned int quantity = BStrToUInt(protein.GetSchema1(QUANTITY_OFF, QUANTITY_LGH));
		++quantity;
		protein.ReplaceSchema1(UIntToBStr(quantity, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
#ifdef LOG_COMPUTATION
		model_time += rate;
		if (type_a) quantity_A.push_back(pair<float, unsigned int>(model_time, quantity + BStrToUInt(protein.GetSchema2(QUANTITY_OFF, QUANTITY_LGH)) + bound));
		else quantity_B.push_back(pair<float, unsigned int>(model_time, quantity + BStrToUInt(protein.GetSchema2(QUANTITY_OFF, QUANTITY_LGH))));
#endif
		return true;
	}
	return false;
}

bool BindUnbindProteins (SCSystem & protein_A, const SCSystem &, SCSystem & protein_B) {
	// Get A and B available quantities
	unsigned int quantity_A = BStrToUInt(protein_A.GetSchema1(QUANTITY_OFF, QUANTITY_LGH)),
				 quantity_B = BStrToUInt(protein_B.GetSchema1(QUANTITY_OFF, QUANTITY_LGH));
	// If some proteins can be bound, bind pairs with the given probability
	unsigned int //max_quantity = MIN(quantity_A, quantity_B),
				 bound_quantity = 0;
	bool is_bound;
	for (unsigned int i=0; i<quantity_A; ++i) {
		is_bound = false;
		for (unsigned int j=bound_quantity; (!is_bound) && (j<quantity_B); ++j) {
			if (frand() <= AB_BIND_RATE) {
				++bound_quantity;
				is_bound = true;
			}
		}
	}
/*
	for (unsigned int i=0; i<max_quantity; ++i)
		if (frand() <= AB_BIND_RATE) ++bound_quantity;
		*/
	unsigned int bound_quantity_A = BStrToUInt(protein_A.GetSchema2(QUANTITY_OFF, QUANTITY_LGH)),
				 bound_quantity_B = BStrToUInt(protein_B.GetSchema2(QUANTITY_OFF, QUANTITY_LGH));
	if (bound_quantity > 0) {
		bound_quantity_A += bound_quantity;
		bound_quantity_B += bound_quantity;
		protein_A.ReplaceSchema2(UIntToBStr(bound_quantity_A, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
		protein_B.ReplaceSchema2(UIntToBStr(bound_quantity_B, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
		quantity_A -= bound_quantity;
		quantity_B -= bound_quantity;
		protein_A.ReplaceSchema1(UIntToBStr(quantity_A, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
		protein_B.ReplaceSchema1(UIntToBStr(quantity_B, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
		protein_A.Grab(protein_B);
		protein_B.Grab(protein_A);
#ifdef LOG_COMPUTATION
		//++nb_computations;
		model_time += AB_BIND_RATE * bound_quantity;
#endif
		return true;
	}
	else if (((bound_quantity_A == 0) || (bound_quantity_B == 0)) && protein_A.HasSubSystem(protein_B) && protein_B.HasSubSystem(protein_A)) {
		protein_A.Release(protein_B);
		protein_B.Release(protein_A);
//#ifdef LOG_COMPUTATION
		//++nb_computations;
//#endif
		return true;
	}
	return false;
}

bool DeleteProteins (SCSystem & protein, const SCSystem &, SCSystem &) {
	// Get protein type
	bool type_a = (protein.GetKernel(K_A_OFF, K_A_LGH) == "0");
	
	// Get degradation rate
	float rate = (type_a) ? A_DEGRADE_RATE : B_DEGRADE_RATE;
	// Get quantity
	unsigned int quantity = BStrToUInt(protein.GetSchema1(QUANTITY_OFF, QUANTITY_LGH));
	// Degrade proteins with the given probability
	unsigned int degraded_quantity = 0;
	for (unsigned int i=0; i<quantity; ++i)
		if (frand() <= rate) ++degraded_quantity;
	if (degraded_quantity > 0) {
		quantity -= degraded_quantity;
		protein.ReplaceSchema1(UIntToBStr(quantity, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
	}
	
	// Get bound quantity
	unsigned int bound_quantity = BStrToUInt(protein.GetSchema2(QUANTITY_OFF, QUANTITY_LGH));
	// Degrade bound proteins with the given probability
	unsigned int unbound_quantity = 0;
	for (unsigned int i=0; i<bound_quantity; ++i)
		if (frand() <= AB_DEGRADE_RATE) ++unbound_quantity;
	if (unbound_quantity > 0) {
		bound_quantity -= unbound_quantity;
		protein.ReplaceSchema2(UIntToBStr(bound_quantity, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
	}

	if ((degraded_quantity > 0) || (unbound_quantity > 0)) {
#ifdef LOG_COMPUTATION
		//++nb_computations;
		model_time += (rate * degraded_quantity) + (AB_DEGRADE_RATE * unbound_quantity);
		if (type_a)	quantity_A.push_back(pair<float, unsigned int>(model_time, bound_quantity + quantity + bound));
		else quantity_B.push_back(pair<float, unsigned int>(model_time, bound_quantity + quantity));
#endif
		return true;
	}
	return false;
}

bool BindProteinToGene (SCSystem & protein, const SCSystem &, SCSystem & gene) {
	// Get protein quantity
	unsigned int quantity = BStrToUInt(protein.GetSchema1(QUANTITY_OFF, QUANTITY_LGH));
	// Bind with the given probability
	bool bind = false;
	for (unsigned int i=0; (i<quantity) && (!bind); ++i)
		if (frand() <= bA_BIND_RATE) bind = true;
	if (bind) {
		gene.ReplaceSchema1("1", INHIBITED_OFF, INHIBITED_LGH);
		--quantity;
		protein.ReplaceSchema1(UIntToBStr(quantity, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
		protein.Grab(gene);
		gene.Grab(protein);
#ifdef LOG_COMPUTATION
		//++nb_computations;
		model_time += bA_BIND_RATE;
		bound = 1;
#endif
		return true;
	}
	return false;
}

bool UnbindProteinFromGene (SCSystem & protein, const SCSystem &, SCSystem & gene) {
	if (frand() <= bA_UNBIND_RATE) {
		// Get protein quantity
		unsigned int quantity = BStrToUInt(protein.GetSchema1(QUANTITY_OFF, QUANTITY_LGH));
		// Remove inhibition
		gene.ReplaceSchema1("0", INHIBITED_OFF, INHIBITED_LGH);
		++quantity;
		protein.ReplaceSchema1(UIntToBStr(quantity, QUANTITY_LGH), QUANTITY_OFF, QUANTITY_LGH);
		protein.Release(gene);
		gene.Release(protein);
#ifdef LOG_COMPUTATION
		//++nb_computations;
		model_time += bA_UNBIND_RATE;
		bound = 0;
#endif
		return true;
	}
	return false;
}

////////////////////////////////////////////////////////////////////////////////
